Centro de Investigación en Economía y Finanzas - CIENFI
Workshop 01: Introducción a Datos Espaciales en R
Eduard Fernando Martínez González
eduard-martinez.github.io
eduard-martinez
efmartinez@icesi.edu.co[1.] Imágenes satelitales
1.1 Motivación
Night lights and economic growth
Asian Financial Crisis: Java, Indonesia
Tomado de: Measuring Economic Growth from Outer
Space
Genocide Event: Rwanda
Tomado de: Measuring Economic Growth from Outer
Space
Puede acceder al documento aquí: Measuring Economic Growth from Outer Space
Night lights and Covid-19
Puede acceder al documento aquí: COVID19 y actividad económica: demanda de energía y luminosidad ante la emergencia
Puede acceder al documento aquí: Examining the economic impact of COVID-19 in India through daily electricity consumption and nighttime light intensity
Measuring the size and growth of cities using nighttime light
Puede acceder al documento aquí: Measuring the size and growth of cities using nighttime light
Housing Vacancy Rate (VHR) y Luces Nocturnas
Puede acceder al documento aquí: An estimation of housing vacancy rate using NPP-VIIRS night-time light data and OpenStreetMap data
CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
Urban Land in Miami, FL
Tomado de: CAUSES OF SPRAWL: A PORTRAIT FROM
SPACE
Puede acceder al documento aquí: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
1.2 ¿Datos raster?
¿Qué es un raster?
Tomado de: https://www.neonscience.org
Resolución
Tomado de: https://www.neonscience.org
Bandas de un raster
Tomado de: https://www.neonscience.org
1.3 Fuentes:
1.3.1 Night lights
DMSPL (1992-2013)
Puede acceder a los datos anuales aquí 1992-2103
VIIRS (2012-Currently)
Puede acceder a la web del proyecto aquí: link
Puede acceder a los datos mensuales aquí: link
Puede acceder a los datos diarios aquí: link
Harmonized (1992-2020)
Puede acceder a los datos y al paper aquí: link
1.3.2 Land cover
1.3.3 NASA
Puede acceder a datos de lluvias, temperatura, contaminación, vientos y otros en la página oficial del proyecto GES DISC de la NASA aquí: https://disc.gsfc.nasa.gov.
[2.] Manipular datos raster en R
0.0 Configuración inicial
Para replicar las secciones de esta clase, primero debe descargar el
siguiente proyecto
de R y abrir el archivo clase-05.Rproj.
Instalar/llamar las librerías de la clase
## Instalar/llamar las librerías de la clase
require(pacman)
p_load(tidyverse,rio,skimr,viridis,osmdata,
raster,stars, ## datos raster
ggmap, ## get_stamenmap
spatstat, ## analis de clusters
sf, ## Leer/escribir/manipular datos espaciales
leaflet) ## Visualizaciones dinámicas##
## The downloaded binary packages are in
## /var/folders/gd/6_fhny_d3y3c9l263sjz73dw0000gn/T//Rtmp4tzXcv/downloaded_packages## solucionar conflictos de funciones
select = dplyr::select2.1 Leer un raster
methods(class = "stars")## [1] [ [[<- [<- %in%
## [5] $<- adrop aggregate aperm
## [9] as_tibble as.data.frame as.im as.owin
## [13] as.POSIXct c coerce contour
## [17] cut dim dimnames dimnames<-
## [21] droplevels expand_dimensions filter hist
## [25] image initialize is.na Math
## [29] merge mutate Ops plot
## [33] prcomp predict print pull
## [37] rename replace_na select show
## [41] slice slotsFromS3 split st_apply
## [45] st_area st_as_sf st_as_sfc st_as_stars
## [49] st_bbox st_coordinates st_crop st_crs
## [53] st_crs<- st_dimensions st_dimensions<- st_downsample
## [57] st_extract st_geometry st_geotransform st_geotransform<-
## [61] st_interpolate_aw st_intersects st_join st_mosaic
## [65] st_normalize st_redimension st_rotate st_sample
## [69] st_set_bbox st_transform st_write time
## [73] transmute write_stars
## see '?methods' for accessing help and source code## importar raster de luces: raster
luces_r = raster('input/raster/night_light_202301.tif')
luces_r## class : RasterLayer
## dimensions : 2990, 2919, 8727810 (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667 (x, y)
## extent : -79.01042, -66.84792, 0.002082733, 12.46042 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : night_light_202301.tif
## names : night_light_202301
## values : -0.64, 2208.49 (min, max)## importar raster de luces: stars
luces_s = read_stars("input/raster/night_light_202301.tif")
luces_s## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## night_light_202301.tif 0.07 0.16 0.18 0.1810952 0.2 0.53 97553
## dimension(s):
## from to offset delta refsys point x/y
## x 1 2919 -79.01 0.004167 WGS 84 FALSE [x]
## y 1 2990 12.46 -0.004167 WGS 84 FALSE [y]2.2 Atributos de un raster
## resolucion
0.00416667*110000 ## [1] 458.3337## geometria
st_bbox(luces_s)## xmin ymin xmax ymax
## -79.010415858 0.002082733 -66.847915761 12.460416166st_crs(luces_s)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]st_dimensions(luces_s)## from to offset delta refsys point x/y
## x 1 2919 -79.01 0.004167 WGS 84 FALSE [x]
## y 1 2990 12.46 -0.004167 WGS 84 FALSE [y]## atributos
names(luces_s)## [1] "night_light_202301.tif"names(luces_s) = "date_202301"
## valores del raster
luces_s[[1]] %>% max(na.rm = T)## [1] 2208.49luces_s[[1]] %>% min(na.rm = T)## [1] -0.64luces_s[[1]] %>% as.vector() %>% summary() ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1 0 0 0 0 2208 4031249luces_s[[1]][is.na(luces_s[[1]])==T] %>% head()# Reemplazar NA's## [1] NA NA NA NA NA NAluces_s[[1]][2000:2010,2000:2010]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21luces_s[[1]][2000:2010,2000:2010] %>% table() # Sustraer una parte de la matriz## .
## 0.0599999986588955 0.0700000002980232 0.0799999982118607 0.0900000035762787
## 3 3 4 3
## 0.100000001490116 0.109999999403954 0.119999997317791 0.129999995231628
## 4 10 10 15
## 0.140000000596046 0.150000005960464 0.159999996423721 0.170000001788139
## 8 14 12 7
## 0.180000007152557 0.189999997615814 0.200000002980232 0.209999993443489
## 6 9 3 4
## 0.219999998807907 0.239999994635582 0.25 0.280000001192093
## 2 2 1 1## puedo reproyectar un raster?
st_crs(luces_s)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]luces_new_crs = st_transform(luces_s,crs=4126)
luces_s[[1]][2000:2010,2000:2010] # no se alteran las geometrias## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21luces_new_crs[[1]][2000:2010,2000:2010] # no se alteran las geometrias## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21## plot data
plot(luces_s)
2.3 Filtrar usando la gemoetría
## download boundary
quilla <- getbb(place_name = "Barranquilla, Colombia",
featuretype = "boundary:administrative",
format_out = "sf_polygon") %>% .[1,]
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla)## cliping
l_quilla_1 = st_crop(x = luces_s , y = quilla) # crop luces de Colombia con polygono de Medellin
l_quilla_1[[1]] %>% t()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA NA NA NA NA NA 6.03
## [26,] NA NA NA NA NA NA NA NA NA NA 2.96 3.27
## [27,] NA NA NA NA NA NA NA NA NA NA 2.14 2.26
## [28,] NA NA NA NA NA NA NA NA NA NA NA 1.99
## [29,] NA NA NA NA NA NA NA NA NA NA NA 1.95
## [30,] NA NA NA NA NA NA NA NA NA NA NA NA
## [31,] NA NA NA NA NA NA NA NA NA NA NA 2.07
## [32,] NA NA NA NA NA NA 2.20 3.50 NA NA 2.04 2.17
## [33,] NA NA NA NA NA 1.65 1.94 1.93 1.88 2.04 2.21 2.26
## [34,] NA NA NA NA 1.78 2.46 2.97 2.20 2.08 3.82 5.71 5.30
## [35,] NA NA NA NA NA 2.48 3.33 2.79 3.26 6.79 9.48 9.97
## [36,] NA NA NA NA NA NA 4.50 8.17 12.00 13.49 11.40 9.58
## [37,] NA NA NA NA 17.51 15.10 21.29 22.38 19.02 11.72 7.46 7.91
## [38,] NA 1.84 9.53 15.32 25.18 18.67 24.94 22.62 14.00 4.76 3.73 3.57
## [39,] NA 4.45 11.74 NA NA 8.18 NA 8.98 6.40 3.54 3.17 3.00
## [40,] NA NA NA NA NA NA NA NA 3.81 3.43 2.76 2.53
## [41,] NA NA NA NA NA NA NA NA NA NA NA NA
## [42,] NA NA NA NA NA NA NA NA NA NA NA NA
## [43,] NA NA NA NA NA NA NA NA NA NA NA NA
## [44,] NA NA NA NA NA NA NA NA NA NA NA NA
## [45,] NA NA NA NA NA NA NA NA NA NA NA NA
## [46,] NA NA NA NA NA NA NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA NA NA NA NA NA NA
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA 0.45 NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA 0.64 NA NA NA NA NA NA
## [6,] NA NA NA NA NA 0.56 NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA 0.78 NA NA NA NA NA
## [10,] NA NA NA NA NA NA 0.86 NA NA NA NA NA
## [11,] NA NA NA NA NA NA 0.89 NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA 1.04 NA NA NA NA
## [13,] NA NA NA NA NA NA 1.21 1.21 1.34 NA NA NA
## [14,] NA NA NA NA NA 1.20 NA NA 2.18 NA NA NA
## [15,] NA NA NA NA NA NA NA NA NA 3.12 NA NA
## [16,] NA 1.60 NA NA NA NA NA NA NA 5.09 14.52 30.40
## [17,] 3.30 6.73 NA NA NA NA NA NA NA NA 32.05 72.30
## [18,] 13.80 16.35 NA NA NA NA NA NA NA 42.66 NA 94.43
## [19,] 32.75 40.53 34.68 NA 12.67 6.47 NA NA NA NA 95.34 76.51
## [20,] 25.76 44.97 46.67 NA NA NA NA NA NA NA 62.60 51.57
## [21,] NA NA NA NA NA NA NA NA NA 15.24 21.01 27.16
## [22,] NA NA NA NA NA NA NA NA NA 40.83 47.05 64.18
## [23,] NA NA 40.23 NA 27.55 NA NA NA 75.98 71.70 85.57 105.79
## [24,] NA 13.75 14.97 15.46 11.35 18.27 29.89 53.61 82.99 93.33 101.92 111.54
## [25,] 5.50 4.56 4.78 4.01 4.50 5.21 14.28 42.38 76.66 85.32 91.35 99.35
## [26,] 3.04 2.84 3.02 3.31 3.98 4.93 14.24 34.13 62.03 71.77 79.13 87.62
## [27,] 2.28 2.48 2.90 3.18 4.69 8.84 23.25 28.93 35.22 45.27 62.79 76.90
## [28,] 2.11 2.35 2.58 3.02 11.18 32.95 71.85 63.93 44.63 33.76 42.14 58.06
## [29,] 2.07 2.20 2.53 3.29 12.92 43.91 82.88 75.04 43.40 22.66 35.22 49.98
## [30,] 2.11 2.26 2.52 2.90 6.60 25.38 28.53 29.68 28.91 42.55 45.05 52.81
## [31,] 2.25 2.42 2.76 2.81 4.35 11.00 16.23 17.50 42.77 62.58 64.82 55.21
## [32,] 2.31 2.49 3.30 4.51 8.08 20.60 36.12 36.38 49.91 70.47 69.25 54.57
## [33,] 2.48 2.82 4.23 7.03 13.51 24.19 39.68 49.92 55.94 69.37 62.20 45.75
## [34,] 5.34 6.25 7.87 10.68 21.66 23.54 28.73 51.42 59.80 58.90 48.16 43.53
## [35,] 12.41 17.60 15.63 22.35 33.51 28.47 29.17 47.31 63.17 52.35 46.63 45.12
## [36,] 15.48 23.54 34.28 43.02 56.55 37.83 52.31 66.26 65.72 53.72 47.85 49.67
## [37,] 11.54 18.39 37.05 68.12 61.43 46.70 58.65 70.88 66.73 56.34 52.74 53.29
## [38,] 4.42 7.21 17.60 36.16 54.00 38.88 30.27 26.63 40.38 57.27 55.93 56.22
## [39,] 3.14 3.63 9.59 16.15 18.85 13.85 10.31 12.22 26.62 54.40 56.93 57.73
## [40,] 2.69 3.13 4.20 7.02 6.92 9.13 12.85 18.41 26.70 60.24 69.57 64.83
## [41,] 2.67 5.43 8.80 8.01 9.52 17.38 25.22 37.87 40.71 54.39 69.89 67.61
## [42,] NA NA NA 39.28 NA NA 31.34 27.85 32.55 47.97 61.82 63.79
## [43,] NA NA NA NA NA NA 17.96 7.92 10.87 25.88 43.63 49.79
## [44,] NA NA NA NA NA NA NA 6.99 8.95 30.74 47.43 54.56
## [45,] NA NA NA NA NA NA NA NA 10.48 28.05 47.61 58.19
## [46,] NA NA NA NA NA NA NA NA 6.33 13.28 27.71 47.00
## [47,] NA NA NA NA NA NA NA NA NA NA 12.52 35.70
## [48,] NA NA NA NA NA NA NA NA NA NA NA NA
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA NA NA NA NA NA
## [17,] 44.35 NA NA NA NA NA NA NA NA NA NA
## [18,] 92.54 72.77 38.77 NA NA NA NA NA NA NA NA
## [19,] 56.01 59.52 47.77 39.21 NA NA NA NA NA NA NA
## [20,] 28.10 32.52 63.07 76.49 85.56 53.18 NA NA NA NA NA
## [21,] 46.55 50.28 71.32 83.99 85.30 63.33 52.52 NA NA NA NA
## [22,] 76.44 77.28 78.88 65.12 60.87 54.69 53.39 43.20 NA NA NA
## [23,] 99.95 84.11 79.16 66.23 53.39 47.17 49.70 51.41 35.46 22.53 NA
## [24,] 102.47 76.40 81.48 79.32 66.51 56.38 59.51 62.07 46.97 25.86 10.23
## [25,] 96.30 91.31 86.57 86.37 84.80 69.04 69.33 66.04 51.06 27.09 19.61
## [26,] 95.27 100.51 87.12 86.48 85.51 77.50 65.75 62.99 51.63 51.20 36.54
## [27,] 89.58 96.36 92.16 90.84 87.09 82.91 75.94 66.18 67.95 74.07 70.44
## [28,] 68.38 83.42 87.67 95.82 90.66 76.26 78.83 78.96 75.55 81.83 70.24
## [29,] 60.20 68.46 79.86 100.32 91.54 73.72 75.44 93.24 92.60 77.20 65.92
## [30,] 58.82 61.30 65.44 69.87 71.60 71.90 75.94 88.63 98.86 86.89 71.64
## [31,] 55.67 55.18 55.81 61.82 66.07 71.64 73.51 78.55 92.67 85.76 82.93
## [32,] 48.39 55.60 57.99 59.12 65.07 69.80 67.73 73.24 83.59 81.23 83.38
## [33,] 43.38 48.67 58.43 58.08 60.75 62.08 58.91 66.09 72.25 70.98 76.85
## [34,] 42.76 48.42 56.79 54.78 54.65 60.34 60.78 67.11 64.15 64.04 64.97
## [35,] 50.09 56.28 57.06 53.92 57.21 72.60 79.05 76.00 68.02 64.23 63.44
## [36,] 57.15 60.50 53.60 50.88 54.49 66.49 70.52 66.77 71.46 68.73 68.48
## [37,] 67.56 65.18 53.31 55.93 60.69 63.25 60.49 60.49 71.85 89.37 73.30
## [38,] 63.81 57.03 52.76 59.40 61.73 75.28 80.32 70.60 82.59 84.18 70.54
## [39,] 57.26 51.44 49.65 50.65 51.92 69.87 86.02 78.78 92.53 87.69 74.86
## [40,] 56.72 47.53 44.60 44.62 49.16 59.19 73.68 80.36 92.06 99.48 96.92
## [41,] 56.60 48.08 45.28 47.42 49.30 67.65 71.78 69.61 84.80 96.47 111.94
## [42,] 62.72 59.50 49.69 52.51 58.76 80.24 76.63 67.45 72.62 79.86 97.51
## [43,] 63.71 60.43 52.18 59.41 72.81 93.18 94.09 87.39 NA NA NA
## [44,] 61.69 58.86 60.99 69.98 79.17 96.70 97.20 NA NA NA NA
## [45,] 64.35 67.04 71.71 82.92 91.30 119.18 NA NA NA NA NA
## [46,] 57.68 NA NA NA NA NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA NA NA NA NA NA
## [,36] [,37] [,38] [,39] [,40] [,41]
## [1,] NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA
## [25,] 12.84 NA NA NA NA NA
## [26,] 31.83 16.36 NA NA NA NA
## [27,] 46.73 21.85 NA NA NA NA
## [28,] 36.41 18.83 11.97 NA NA NA
## [29,] 41.21 23.91 22.20 NA NA NA
## [30,] 60.35 49.21 42.73 29.10 NA NA
## [31,] 82.86 69.00 61.91 51.15 NA NA
## [32,] 95.06 82.41 64.15 52.33 45.85 NA
## [33,] 87.11 81.52 44.85 37.96 39.19 NA
## [34,] 68.56 65.19 51.40 72.00 66.22 NA
## [35,] 70.97 66.07 56.48 89.83 82.82 NA
## [36,] 75.41 76.00 55.55 51.73 47.74 NA
## [37,] 72.00 73.48 63.38 50.81 37.12 NA
## [38,] 71.74 76.27 80.48 83.79 80.64 NA
## [39,] 83.80 94.17 93.37 85.02 82.15 NA
## [40,] 87.09 91.57 88.35 76.77 37.66 21.85
## [41,] NA NA NA NA NA NA
## [42,] NA NA NA NA NA NA
## [43,] NA NA NA NA NA NA
## [44,] NA NA NA NA NA NA
## [45,] NA NA NA NA NA NA
## [46,] NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA## Plot data
ggplot() + geom_stars(data=l_quilla_1 , aes(y=y,x=x,fill=date_202301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 
[3.] Aplicación: Actividad económica en Barranquilla
3.1 Importar datos 201301
## load data
l_quilla_0 = read_stars("input/raster/night_light_201301.tif") %>% st_crop(quilla)
names(l_quilla_0) = "date_201301"
ggplot() + geom_stars(data=l_quilla_0 , aes(y=y,x=x,fill=date_201301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 
## stack rasters: adicionar bandas
l_quilla = c(l_quilla_0,l_quilla_1)
l_quilla## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## date_201301 0.15 6.1175 48.555 42.78252 69.3050 135.63 1270
## date_202301 0.45 15.2600 51.175 46.05638 69.8525 119.18 1270
## dimension(s):
## from to offset delta refsys point x/y
## x 982 1022 -79.01 0.004167 WGS 84 FALSE [x]
## y 325 372 12.46 -0.004167 WGS 84 FALSE [y]3.2 Raster a datos vectoriales
## st_as_sf
puntos_quilla = st_as_sf(x = l_quilla, as_points = T, na.rm = T) # raster to sf (points)
puntos_quilla## Simple feature collection with 698 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.92083 ymin: 10.9125 xmax: -74.75417 ymax: 11.10833
## Geodetic CRS: WGS 84
## First 10 features:
## date_201301 date_202301 geometry
## 1 0.15 0.45 POINT (-74.85417 11.10417)
## 2 0.26 0.64 POINT (-74.85 11.09167)
## 3 0.27 0.56 POINT (-74.85 11.0875)
## 4 0.32 0.78 POINT (-74.84583 11.075)
## 5 0.42 0.86 POINT (-74.84583 11.07083)
## 6 0.39 0.89 POINT (-74.84583 11.06667)
## 7 0.92 1.04 POINT (-74.84167 11.0625)
## 8 0.91 1.21 POINT (-74.84583 11.05833)
## 9 1.25 1.21 POINT (-74.84167 11.05833)
## 10 1.13 1.34 POINT (-74.8375 11.05833)poly_quilla = st_as_sf(x = l_quilla, as_points = F, na.rm = T) # raster to sf (polygons)
poly_quilla## Simple feature collection with 698 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.92292 ymin: 10.91042 xmax: -74.75208 ymax: 11.11042
## Geodetic CRS: WGS 84
## First 10 features:
## date_201301 date_202301 geometry
## 1 0.15 0.45 POLYGON ((-74.85625 11.1062...
## 2 0.26 0.64 POLYGON ((-74.85208 11.0937...
## 3 0.27 0.56 POLYGON ((-74.85208 11.0895...
## 4 0.32 0.78 POLYGON ((-74.84792 11.0770...
## 5 0.42 0.86 POLYGON ((-74.84792 11.0729...
## 6 0.39 0.89 POLYGON ((-74.84792 11.0687...
## 7 0.92 1.04 POLYGON ((-74.84375 11.0645...
## 8 0.91 1.21 POLYGON ((-74.84792 11.0604...
## 9 1.25 1.21 POLYGON ((-74.84375 11.0604...
## 10 1.13 1.34 POLYGON ((-74.83958 11.0604...## Plot data
ggplot() +
geom_sf(data = puntos_quilla , aes(col=date_202301)) +
scale_color_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test() 
ggplot() +
geom_sf(data = poly_quilla , aes(fill=date_202301),col=NA) +
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test()
3.3 Asignar el valor del pixel a un lugar
## get polygon
puerto <- getbb(place_name = "Puerto de Barranquilla - Sociedad Portuaria, Colombia",
featuretype = "barrier:wall",
format_out = "sf_polygon")
leaflet() %>%
addTiles() %>%
addPolygons(data=puerto)## asignar luminosidad
l_puerto = st_join(puerto,poly_quilla)
## variacion promedio
df = l_puerto
st_geometry(df) = NULL
df %>%
summarise(pre=mean(date_201301,na.rm=T) ,
post=mean(date_202301,na.rm=T)) %>%
mutate(ratio=post/pre-1)## pre post ratio
## 1 52.647 59.537 0.1308717[4.] Aplicación: Cambio de cobertura de suelo?
## read raster
land_cover = c(read_stars("input/raster/discreta_2015.tif"),
read_stars("input/raster/discreta_2019.tif"))
0.000992063*111000 # resolucion## [1] 110.119## read polygon
quilla_buffer <- st_buffer(x = quilla,dist = 2000)
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_buffer)quilla_border <- st_difference(x = quilla_buffer , y = quilla)
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_border)## cliping raster
quilla_cover = land_cover %>% st_crop(quilla_buffer)
plot(quilla_cover) ## view data
## cliping raster
border = land_cover %>% st_crop(quilla_border)
plot(border) ## view data
## raster to sf
border_sf = st_as_sf(x=border, as_points = F, na.rm = T)
## cuantos pixeles se urbanizaron?
border_sf = border_sf %>%
rename(c_2015=discreta_2015.tif,c_2019=discreta_2019.tif)
border_sf = border_sf %>%
mutate(build = ifelse(c_2015!=50 & c_2019==50,1,0))
table(border_sf$build)##
## 0 1
## 13990 46[5.] Aplicación: Clusteres espaciales
## restaurantes
points <- opq(bbox = getbb("Bogota Colombia")) %>%
add_osm_feature(key = "amenity", value = "restaurant") %>%
osmdata_sf() %>% .$osm_points %>% select(osm_id,name) %>% mutate(restaurant=1)
leaflet() %>% addTiles() %>% addCircles(data=points)## download boundary
city <- getbb(place_name = "Bogota, Colombia",
featuretype = "boundary:administrative",
format_out = "sf_polygon") %>% .[[2]]
leaflet() %>% addTiles() %>% addPolygons(data=city)## afine transformations
points = st_transform(points,3857)
city = st_transform(city,3857)
## sf to ppp
points_p <- as.ppp(X = points["restaurant"])
Window(points_p) <- as.owin(city)
plot(points_p)
#--------------#
## summary
summary(points_p) # las unidades son en metros## Marked planar point pattern: 6882 points
## Average intensity 1.733022e-05 points per square unit
##
## Coordinates are given to 10 decimal places
##
## marks are numeric, of type 'double'
## Summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 7590 395161000.0 0.995000
## polygon 2 276 1524740.0 0.003840
## polygon 3 72 235244.0 0.000592
## polygon 4 45 75844.4 0.000191
## polygon 5 81 112582.0 0.000284
## enclosing rectangle: [-8262524, -8238784] x [498234.9, 538674.6] units
## (23740 x 40440 units)
## Window area = 397110000 square units
## Fraction of frame area: 0.414intensity(points_p)## [1] 1.733022e-05## reescalar
points_p = rescale(X = points_p , s = 100 , unitname = "100-metros")
unitname(points_p)## 100-metros / 100-metrossummary(points_p)## Marked planar point pattern: 6882 points
## Average intensity 0.1733022 points per square 100-metros
##
## Coordinates are given to 12 decimal places
##
## marks are numeric, of type 'double'
## Summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 7590 39516.10000 0.995000
## polygon 2 276 152.47400 0.003840
## polygon 3 72 23.52440 0.000592
## polygon 4 45 7.58444 0.000191
## polygon 5 81 11.25820 0.000284
## enclosing rectangle: [-82625.24, -82387.84] x [4982.349, 5386.746] 100-metros
## (237.4 x 404.4 100-metros)
## Window area = 39711 square 100-metros
## Unit of length: 1 100-metros
## Fraction of frame area: 0.414intensity(points_p)## [1] 0.1733022## estimated standard error for intensity
sqrt(intensity(points_p)/summary(points_p)$window$area)## [1] 0.00208904## plot
plot(density(points_p,5)) ## plot density
contour(density(points_p,bw.ppl(points_p)), axes = FALSE) ## contour density
#--------------#
## quadrants
points_q = quadratcount(points_p, nx = 4, ny = 6)
points_q## tile
## Tile row 1, col 3 Tile row 1, col 4 Tile row 2, col 2 Tile row 2, col 3
## 0 20 100 345
## Tile row 2, col 4 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3
## 739 0 373 861
## Tile row 3, col 4 Tile row 4, col 1 Tile row 4, col 2 Tile row 4, col 3
## 864 53 647 2029
## Tile row 4, col 4 Tile row 5, col 1 Tile row 5, col 2 Tile row 5, col 3
## 263 17 183 377
## Tile row 5, col 4 Tile row 6, col 2 Tile row 6, col 3
## 0 0 11plot(points_q, cex = 2 , col="red")
Referencias
Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]
- Cap. 4: Spatial Data Operations
- Cap. 5: Geometry Operations
- Cap. 6: Reprojecting geographic data
- Cap. 11: Statistical learning
Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]
- Cap. 4: Spatial Data Import and Export
- Cap. 7: Spatial Point Pattern Analysis
- Cap. 8: Interpolation and Geostatistics